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Mergers and ejections of black holes in globular clusters 
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ABSTRACT 

We report on results of fully consistent TV-body simulations of globular cluster models with 
TV = 100 000 members containing neutron stars and black holes. Using the improved 'algo- 
rithmic regularization' method of Hellstrom and Mikkola for compact subsystems, the new 
code NB0DY7 enables for the first time general relativistic coalescence to be achieved for post- 
Newtonian terms and realistic parameters. Following an early stage of mass segregation, a few 
black holes form a small dense core which usually leads to the formation of one dominant bi- 
nary. The subsequent evolution by dynamical shrinkage involves the competing processes of 
ejection and mergers by radiation energy loss. Unless the binary is ejected, long-lived triple 
systems often exhibit Kozai cycles with extremely high inner eccentricity (e > 0.999) which 
may terminate in coalescence at a few Schwarzschild radii. A characteristic feature is that 
ordinary stars as well as black holes and even BH binaries are ejected with high velocities. On 
the basis of the models studied so far, the results suggest a limited growth of a few remaining 
stellar mass black holes in globular clusters. 

Key words: black hole physics - globular clusters: general - methods: numerical 



1 INTRODUCTION 

Recent years have seen many studies relating to the dynamics of 
black holes (BHs) both in galactic and extra-galactic systems. In 
view of the observations of the S stars at the galactic centre it is 
not surprising that most efforts have been directed towards a rel- 
atively massive BH. However, there is also considerable interest 
J* in the effect of BHs in star clusters. It has been argued that be- 
. C-> cause velocity kicks may occur at formation, some BHs are more 
likely to be retained in globular rather than open clusters. We may 
<_i therefore distinguish between problems dealing with several stellar 
mass BHs or one dominant body formed by an accretion process. 
As far as TV-body simulations are concerned, these two types call 
for different methods of solution. Thus in the former we need to 
treat strong interactions of BHs, while in the latter case a num- 
ber of short binary periods require careful attention. A closer in- 
spection of investigations concerned with a relatively massive BH 
shows that Newtonian motions are often adopted in full simulations 
(Brockamp, Baumgardt & Rroupa 201 1). On the other hand, inves- 
tigations involving binary BHs tend to avoid the numerical prob- 
lems associated with small Schwarzschild radii by artificial scaling 
of the masses or magnifying the post-Newtonian effects (Iwasawa, 
Funato & Makino 2006, Berentzen et al. 2009). 

The first full N-body simulation with post-Newtonian (PN) 
terms (Aarseth 2003b) indicated that a black hole binary of inter- 
mediate mass may achieve the general relativistic (GR) coalescence 
condition following eccentricity growth by the Kozai mechanism 
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(Kozai 1962, 2004). This development was made possible by a spe- 
cial integration method for a long-lived massive binary (Mikkola & 
Aarseth 2002, TTL). An alternative method, suitable for a dominant 
single BH, showed several examples of GR coalescence for realis- 
tic parameters (Aarseth 2007). Although expensive, the so-called 
wheel-spoke method generalized from a three-body regularization 
algorithm (Aarseth & Zare 1974; Zare 1974) is accurate and en- 
ables extremely close orbits to be studied. It is notable that in both 
these investigations, Kozai cycles played an important role, with 
maximum eccentricities often reaching large values (e.g. 0. 99990. 

The motivation for the present investigation goes back more 
than 40 years. Thus an early simulation demonstrated that two mas- 
sive bodies (factor of 5) gave up kinetic energy to the other mem- 
bers and reached the centre after a few crossing times (Aarseth 
1971). Since most of the other core particles were expelled, the 
two dominant members invariably formed a binary. Although the 
cluster membership was only TV = 250, the mass segregation still 
operated in a qualitatively similar way. Given a dynamically shrink- 
ing massive binary, the probability of a long-lived triple with suit- 
able inclination is non-negligible, especially since any stellar mass 
suffices for inducing Kozai cycles. 

Further simulations of small globular clusters with N — 10 5 
members containing a significant BH component were made with 
special-purpose GRAPE computers (Mackey et al. 2008). More re- 
cently similar systems have been studied with Graphics Process- 



1 Early TV-body simulations report large eccentricity growth to e > 0.999 
at constant semi-major axis in triple systems which exhibit the hallmarks of 
Kozai cycles (van Albada 1968, Aarseth 1971). 
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ing Units (GPU) which also require special programming (Nitadori 

2009) . The software implementation in the standard NBODY6 code 
based on the fourth-order Hermite block-step neighbour scheme 
has been described in a recent paper (Nitadori & Aarseth 2012). 
Of special interest here is an application to compact subsystems of 
more than two dominant members, all of which may be BHs. In 
particular, we wish to explore possible post-Newtonian dynamics 
in a globular cluster context. Although there have been attempts 
of including the 2.5PN radiation term using the standard NB0DY6 
code (Aarseth 2003a) by others (Banerjee, Baumgardt & Kroupa 

2010) , simplified treatments fall short of realistic requirements and 
hence a fresh approach is called for. 

In this paper, we make use of a new powerful integration pack- 
age that can overcome the main numerical difficulties associated 
with compact subsystems. Here the challenge is to treat stellar mass 
BH binaries of extremely short periods during late stages of in- 
spiral using realistic parameters. The improved algorithmic regu- 
larization chain method (Hellstrom & Mikkola 2010, ARC) traces 
its development from chain regularization (Mikkola & Aarseth 
1993, CHAIN) via the logarithmic Hamiltonian method (Mikkola 
& Tanikawa 1999) and the time-transformed leapfrog code TTL 
to special treatments of post-Newtonian terms (Mikkola & Merritt 
2008), also denoted algorithmic regularization. The resulting code 
called NB0DY7 is able to deal with extreme few-body configura- 
tions up to the final stages of GR coalescence without recourse to 
artificial scaling of the parameters. A key feature of the present in- 
vestigation is that velocity kicks are assigned to neutron stars and 
BHs at formation, such that only about 10 percent of the member- 
ship is retained in both cases. The subsequent evolution by mass 
segregation of the heaviest members gives rise to compact sub- 
systems which sometimes require post-Newtonian treatments due 
to the process of Kozai eccentricity cycles. Consequently, even a 
small population of stellar mass black holes are able to control the 
central region. The challenge is then to employ an efficient method 
for dealing with extremely compact subsystems where PN terms 
may be needed for the dominant two-body motion. 

This paper is organized as follows. We begin by outlining 
the A r -body implementation of algorithmic regularization and post- 
Newtonian terms for a compact subsystem. Section 3 is devoted to 
a general description of mass segregation while results for particu- 
lar models are illustrated in section 4. Finally, some aspects of the 
simulations are discussed in section 5, followed by conclusions. 



2 COMPACT SUBSYSTEMS 

Simple test calculations with a few free-floating stellar mass BHs 
added reveal a characteristic behaviour of mass segregation in ac- 
cordance with expectations. The essential question is concerned 
with the degree of central concentration; namely when does the 
innermost core stop shrinking due to binary activity. Some prelim- 
inary investigations were sufficiently encouraging to proceed with 
new code developments as summarized below. 

The main challenge for dealing with a compact configura- 
tion is to maintain accuracy, albeit at increased cost. This involves 
some kind of regularization to avoid near-singular force terms. The 
adopted method (Mikkola & Merritt 2006, 2008) is based on a so- 
phisticated time transformation that enables even collision orbits 
to be integrated at high accuracy. Briefly, the underlying integrator 
(Bulirsch & Stoer 1966) senses an approach to collision and avoids 
evaluating the force at singular points. Moreover, using a judicious 
choice of coefficients in the time transformation, the problem of 



large mass ratios can also be accommodated. Two analogous forms 
of the time transformations are used for a subsystem of iV c h inter- 
acting particles (Mikkola & Merritt 2008), 

ds = [a(T + B)+/3uj + j]dt = [all + /3fi + j]dt, (1) 

where ds is the new differential time element and a, ft, j are ad- 
justable constants. Furthermore, T is the kinetic energy, U the (pos- 
itive) potential energy, B the binding energy, B = U — T, and Q 
is an optional function of the coordinates. It can be shown that this 
relation provides well-behaved solutions for two-body collisions 
when used in connection with a simple leapfrog algorithm where 
the results are improved by an extrapolation method. Moreover, 
the case of velocity-dependent perturbations (e.g. post-Newtonian 
terms) can be treated explicitly with high accuracy (Hellstrom & 
Mikkola 2010). 

The alternative time transformations are applied to coordi- 
nate and velocity leapfrog integrations, respectively. Thus the Q- 
formulation is related to the TTL method in which cb — v-Vfl 
provides a regular solution for u(t) — fi(t) (Mikkola & Aarseth 
2002). When a range of masses are involved, the choice (a, /3, 7) = 
(1.0,0.001,0) is recommended, where the non-zero value of f3 
yields increased accuracy for any massive bodies. Note that the 
algorithmic chain regularization actually employs the chain data 
structure without solving the corresponding equations of motion. 
This procedure leads to significant reduction of round-off errors 
and therefore plays a key role in the formulation. 

The choice of the subsystem membership is to some extent 
experimental. Thus we must consider the cost of treating a small 
subsystem by an accurate method where a larger size would also 
necessitate a greater number of perturbers. Note that the basic in- 
tegrator performs a large number (~ 100) of function evaluations 
per step and hence also requires many coordinate predictions of the 
perturbers for consistency. Although very generous choices of the 
perturber number have been used before (Harfst et al. 2008), it has 
yet to be demonstrated that this is necessary. A careful investiga- 
tion of the corresponding energy change reveals that the external 
effect is relatively small for dimensionless perturbations beyond 
7 P ert ~ 10~ 9 which is a typical limit for selection (reduced from 
the earlier value 10~ 7 ). Even so, perturber numbers of only 3 or 
4 are typical in the present simulations with small chain member- 
shipfl More problematic are the decisions about when a perturber 
should be included in the membership and vice versa. 

The software package ARC replaces the procedures relating to 
the CHAIN code, with analogous decision-making for the interface 
connection to the standard NBODY6. Consequently, the new code 
is called NBODY7. Although ARC contains post-Newtonian terms 
up to 3.5PN, an alternative formulation up to 3PN (Blanchet & Iyer 
2003, Mora & Will 2004) has been retained (Aarseth 2007). This 
has the advantage that relativistic expressions for the semi-major 
axis and eccentricity are readily available for decision-making and 
data analysis. 

We now review some relevant aspects relating to the imple- 
mentation of PN terms. The equation of dominant two-body motion 
with G = 1 is written as 

where M = mi + m-2 and A and B represent the post-Newtonian 



2 This behaviour is characteristic for a dominant binary. 

© 2011 RAS, MNRAS 000,[T_p?? 



Black Holes in Globular Clusters 3 



terms for coordinates and velocities, r and v, respectively. The per- 
turbing force then takes the form 



GR = 



mitn2 



/2 



(B 1 + ^+ B5/2 



r3 



+ c 4 >r 



(3) 



In iV-body applications, the scaled speed of light is formally given 
by c = 3x 10 5 /V* , with V* the velocity unit (typically 16 km s _1 
here). However, much smaller values are often used by other work- 
ers for easing numerical problems. The expansion coefficients in 
A, B to sixth order in 1/c are obtained from the current formula- 
tion. Although the spin term is of order c~ 3 , it has only been used in 
tests here for the most dominant body (Gopakumar, private commu- 
nication). One reason for this neglect is that the spin effect depends 
on the uncertain initial magnitude as well as direction. However, 
only the radiation term gives rise to energy dissipation. 

The time-scale for gravitational radiation (in iV-body units) is 
given by the classical expression (Peters 1964) 



TGR 



(1 



2\7/2 



(4) 



64X(1 + X)mf ff (e) 

where X = 7712/7711 with m 1 > 7112 and g(e) a known function 
of eccentricity (~ 4.5 for large e). In the relativistic regime, the 
binding energy per unit mass is determined from 



. ei £2 £3 



(5) 



(6) 



where eo is the Newtonian value. This yields the semi-major axis 
- -M. 

Likewise, the eccentricity is obtained via the angular momentum 
expansion 



h , h 



by 



J = J (l + ^ + 



~ (1 Ma' 



(7) 



(8) 



The coefficients ei, ea, £3, /i, fa are also taken from the literature 
cited above. 

Another useful quantity is the Lenz vector 

v x r x v r 



AI 



(9) 



which can be used to determine the Einstein (1915) pericentre shift 
(in radians per orbit) 

67rAf 



Ait; = 



c 2 a(l 



(10) 



An effort has been made to maintain a high level of energy 
conservation. For this purpose we split the contributions to subsys- 
tem energy changes into two parts. Thus the integrated effect of 
particle perturbations is used to monitor the instantaneous internal 
energy while additional contributions arising from the PN terms 
are treated separately. Although only the gravitational radiation is 
dissipative, the other terms depend on the orbital parameters and 
therefore need to be included for conservation purposes. Conse- 
quently the sum of the internal energy and the accumulated rela- 
tivistic energy change is added to the standard iV-body value based 
on considering the subsystem as a point-mass body. 

The question of when to include PN terms is a delicate one. In 



view of the increasing cost of the higher orders, we have devised a 
scheme based on efficiency which has proved itself (Aarseth 2007). 
The basic idea is to delay the PN stage until the radiation time-scale 
falls below a specified value; e.g. tgr, < 500 or about 30 Myr. 
At this stage, the first-order precession terms of comparable cost 
are also activated, especially because a high precession rate may 
affect any eccentricity growth. The second and third-order terms 
are then included once the time-scale shrinks to 50 and 1 iV-body 
units, respectively. Inspection of actual examples show that once or- 
bital shrinkage begins, a significant acceleration often takes place, 
boosted by favourable Kozai cycles. 

As an additional safeguard, the next order is also activated if 
the Kozai period (Kiseleva, Eggleton & Mikkola 1998), Tkozai, 
falls below a certain value. This time-scale can be evaluated for 
long-lived triple configurations which are often present. The ulti- 
mate aim of post-Newtonian iV-body simulations is to see whether 
coalescence can occur for realistic parameters; in this case stellar 
mass BHs in globular clusters. Although even smaller two-body 
separations can be reached with high accuracy, we now define coa- 
lescence at four Schwarzschild radii by 



Rco&l 



8M 
^2 • 



(11) 



As can be seen from equation ©, the corresponding time-scale 
becomes extremely small during the final stage^j. 

Termination by coalescence may also be defined at an earlier 
stage if conditions are favourable. We distinguish between three 
different cases; (i) N c ^ = 2 with no perturbers and small tgr (e.g. 
tgr < 1), (ii) a(l — e) < i? CO ai during the second or third PN 
stage, and (iii) N c ^ = 3, Tkozai > 25, likewise for small tgr, and 
large outer pericentre (factor of 100). 

At later times it frequently happens that a dominant binary 
may not have any perturbers and the GR time-scale is large, in 
which case it is treated in the usual unperturbed two-body approx- 
imation; i.e. with the internal motion frozen. A new subsystem is 
then re-initialized (at the same orbital phase) once a binary per- 
turber is identified, based on several conservative procedures in- 
volving relative distances and velocities (Aarseth 2003a, p,190)Q. 

The ARC system is advanced in a similar way to the standard 
CHAIN code. Thus each block-step time interval is treated sequen- 
tially before control is returned to the main code which also deals 
with regularized binaries in a similar manner. A full discussion of 
the GPU implementation in NBODY6 has been presented elsewhere 
(Nitadori & Aarseth 2012). 



3 MASS SEGREGATION 

We have performed a number of iV-body simulations in order to 
improve the statistical results of rare events. Equilibrium Plummer 
models are used with N = 1 x 10 5 particles and an IMF mostly 
in the range 50 - 0.1 M (Kroupa, Tout & Gilmore 1993) ex- 
tended into the brown dwarf regime. We adopt standard iV-body 
units with G = 1 where the total mass and energy are scaled 
to 1 and E = —1/4, respectively, with mean square velocity of 
1/2. Physical units are then readily obtained once the length scale 
(R* in pc) and mean mass (m = 1/N in Mq) are specified. For 
most of the present models R* — 1 pc and rh = 0.6 Mq, with 



3 Typical parameters are a = 1 X 10~ 10 , c = 2 X 10 4 , mi = 3 X 10 4 . 

4 Following initialization, the relative perturbation is typically < 10 — 7 . 
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2000 4000 6000 6000 10000 12000 14000 

Time 

Figure 1. Mean central distance for light masses (upper curve), heavy 
masses (middle curve) and BHs (lower curve), Model BH7. Because 
of mass-loss evolution, the maximum stellar mass has been reduced to 
3.1 Mq at the end (750 Myr), with the mid-point and lower values at 
0.61 Mq and 0.10 Mq , respectively. A small number of BHs are included 
in the heavy population but contribute reduced values. 



a scaled half-mass radius ~ 0.8 and time-scale conversion fac- 
tor T* = 0.06 Myr. Likewise, the typical velocity scale unit is 
V* ~ 16 km s -1 . Since the speed of light is traditionally taken as 
a free parameter in such simulations, we have adopted c = 18 000 
for the final models, close to the actual value. 

Before presenting results, it may be of interest to consider 
the justification for the post-Newtonian implementation in the con- 
text of the present parameters. Let us assume an energetic binary, 
formed by dynamical means with semi-major axis ahard — 2/N = 
2 x 10~ J and component masses 20 Mq. The energy of the mas- 
sive binary would then represent 1 percent of the total energy which 
can be achieved in this type of simulation. Moreover, from previ- 
ous experience (Aarseth 2007), a maximum eccentricity of 0.999 
is often seen. According to equation {4}, the corresponding time- 
scale would be 35 A^-body units or 2 Myr which is practically fea- 
sible, although it would require a very large number of binary or- 
bits. Moreover, starting the PN sequence at such high eccentricity 
usually implies the presence of a Kozai cycle which would tend to 
counteract the decay and hence assist in reducing the actual time 
interval further. 

The full A^-body simulations employ the NBODY6 synthetic 
stellar evolution package (Hurley, Pols & Tout 2000, 2002) which 
provides information about stellar mass, radius and type as a func- 
tion of time. Of particular relevance here are the final BH masses 
(Eldridge & Tout 2004). Thus for an intermediate metallicity of 
Z = 0.001 each model typically yields about 800 neutron stars 
and 140 BHs. The latter are in the range 3.0 to 19 Mq with an 
average mass of 8 Mq . In the absence of a consensus on BH and 
neutron star velocity kicks at formation, we adopt a conservative 
approach. Thus the kick distribution is chosen from a Maxwellian 
with dispersion 2V* ~ 32 kms -1 and applied to both popula- 
tions, which gives a retention of about 10 percent. As usual with 
tidal fields, escaping stars are removed at twice the tidal radius, 
where r t id c — 56 pc. 

The increased velocity dispersion for BHs due to velocity 
kicks is evident during the time interval of 4 — 10 Myr, followed 
by neutron star formation up to 60 Myr. A further stage of escape 
reduces the populations to about 14 and 40, respectively, with an 
average BH mass of 9 Mq at age 100 Myr. By this time, the mass 



Model BH8 M1 
Model BH8 M2 




Figure 2. Mean central distance for light masses (upper curve), heavy 
masses (middle curve) and BHs (lower curve), Model BH8. The mass 
ranges are similar to Fig. 1 . 

segregation of BHs is well developed. A typical behaviour is exhib- 
ited in Fig. 1 which shows the mean central distance of two popula- 
tions, each containing 50 percent of the total mass for Model BH7. 
The third plot gives the similar information for the BHs. Since the 
BH membership is relatively small, the geometrical mean separa- 
tion is used here. This definition is more representative when the 
result may be affected by a few outliers, including escapers. Un- 
less specified otherwise, all quantities plotted are in A^-body units, 
where the scaling factor for time in Myr is 0.06. 

We note the characteristic behaviour that the mean central 
distance for the heavy members also increases with tim^f). The 
small BH population exhibits significant contraction subject to spo- 
radic expansion. The rapid expansion of the core from small values 
(~ 0.05) near the times 3200 and 6300 is connected with signif- 
icant binary activity. In fact, the first BH binary is created at the 
former time (a ~ 5 x 10~ 4 ). Further shrinkage of the expanded 
core then takes place until a similar minimum is reached where- 
upon an existing binary experiences exchange with fast ejection 
(« j ~ 48 km s~ x ) of a massive (11 Mq) component. It is worth 
noting that already at t ~ 3000 (or 180 Myr) only 14 BHs re- 
main. A third notable event occurs near t ~ 9600 when the binary 
is ejected in a slingshot interaction from a bound triple. The rela- 
tive centre of mass (cm.) velocities of the ejected members, 20 and 
63 km s~\ mirror closely the mass ratio of 3.2. Thus we see four 
distinct stages of core collapse, with only six relatively light BHs 
remaining when the calculations were terminated. 

The GR radiation time-scales of ejected BH binaries are of 
potential interest for future detectors. In the absence of a systematic 
study, we note that nearly all the values predicted by equation $4$ 
are well below 10 Gyr with 10-100 Myr being typical, in qualitative 
agreement with provisional findings (Banerjee et al. 2010). 

More pronounced evidence of oscillatory core behaviour 
can be seen in Fig. 2 which displays the same quantities for 
Model BH8. Here a genuine GR coalescence occurs near t ~ 5060 
but this does not induce core expansion. However, significant bi- 
nary activity occurs around t ~ 5950 which results in the fast 
ejection of a massive BH and subsequent expansion. The next no- 
ticeable expansion is at t ~ 8970 when a massive triple escaped 
by the slingshot mechanism. Again the large ejection velocities of 
29 and 121 km s _1 with respective masses 45 and 11 Mq confirm 

5 The contrast is enhanced further using the geometric mean separation. 
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Figure 3. Geometric mean central distance for neutron stars (upper curve) 
and black holes (lower curve), Model BH8. 



Figure 4. Ejection velocities in km s 1 from the central BH binary, Model 
BH7. 



nicely the momentum conservation. It should be emphasized that 
the slingshot events above were purely Newtonian, with the final 
transition from a ~ 8 x 10~ 6 to 4 x 1CT 6 . 

A comparison between neutron stars and BHs for the last 
model is shown in Fig. 3, where the mean distance for the neu- 
tron stars is also evaluated using the geometric mean. We note the 
steady trend for the final 41 neutron stars to segregate towards the 
centre, with only six remaining BHs of mean mass 4 Mq. Thus 
on a slightly longer time-scale the core would be dominated by the 
neutron stars, nearly all of mass 1.4 Mq . This difficult and long cal- 
culation consumed nearly 300 hours but even so, the accumulated 
energy error of 7 x 10 is highly satisfactory for such a demand- 
ing integration. Finally, the general expansion of the whole cluster 
amounts to a factor of 3 as measured by the half-mass radius, in 
agreement with a standard calculation (Gieles et al. 2010, eq. 6), 
while the actual binding energy is reduced by 8. Some of the early 
expansion of the massive component is undoubtedly due to the ef- 
fects of mass loss from evolving stars. It may also be noted that 
the mass segregation and formation of a compact subsystem due to 
dynamical friction is in accordance with earlier theoretical analysis 
(Spitzer 1969) although a direct comparison would be difficult. 



4 MODELS 

In this paper, we are mainly interested in the behaviour of black 
holes. Before looking at detailed results, we summarize some char- 
acteristic stages in a selected model displayed below, where i de- 
notes the inclination. Here the key word "Subsystem" means that 
at least two BHs are sufficiently compact for the new treatment. 

• BH binary, t = 465, N bh = 9, a = 7 x 10" 5 

• Subsystem, t = 446, a = 3 x 10~ 5 , e = 0.75 

• Shrinkage, t = 523, a = 1.8 x 10~ 5 , e = 0.65 

• Kozai cycles, t = 552, i = 98, Tkozai = 0.1 Myr 

• Eccentricity, t = 553, a = 8 x 10~ 6 , e = 0.9996 

• Coalescence, t = 554, a = 1 x 10~ n , e = 0.03 

• Eccentricity, t = 564, a = 3 x 10~ 4 , e = 0.9999 

• Shrinkage, t = 638, a = 3 x 10" 5 , e = 0.67 

• Subsystem, t = 718, a = 2 x 10" 5 , e = 0.999 

• Coalescence, t = 720, a = 4 x 10~ 8 , tgh, = 1 

A few key features in Model BH12 are displayed at different 
times in Myr. This model is somewhat unusual in that there are 



two coalescence events, with the most massive BH formed by suc- 
cessive mergers, first by combining 11.9 and 13.4 Mq, followed 
by the accretion of 15.6 Mq. In spite of continuing the calcula- 
tion another 700 Myr, the massive BH remained in the system be- 
cause during this stage the total mass of the other BH members 
only amounted to 19 Mq. We note that each coalescence was pre- 
ceded by high eccentricity which resulted in time-consuming post- 
Newtonian calculations until the condition Jilt was reached or the 
time-scale became sufficiently short. 

The central concentration of BHs eventually leads to strong in- 
teractions. Particularly energetic ejections occur after the formation 
of a dominant binary and at some stage even the binary may escape 
due to the recoil effect. All stars with high ejection velocities in 
Model BH7 are shown in Fig. 4. Here the velocity is evaluated with 
respect to the subsystem cm., usually a binary, taking the lower 
limit of twice the current average velocity which typically ensures 
escape. In this model, the binary itself was ejected by recoil at time 
9600 with velocity i; c j = 20 km s _1 and actual velocity at escape 
"esc = 10 km s _1 . There were four other BHs with terminal es- 
cape velocity above 20 km s _1 but no dynamically ejected neutron 
stars above half this valu^fl In spite of these strong interactions, 
the numerical accuracy was maintained well with the accumulated 
energy error only amounting to 2 x 10 _J . It is also worth noting 
that the average BH mass declines to about 6 Mq which is a direct 
consequence of mass segregation where the most massive BHs are 
ejected first. 

Although energetic, the ejection of massive binaries invariably 
occurs using the Newtonian formulation where the semi-major axis 
is still within a modest factor of Ohard- Consequently, favourable 
configurations are initiated through the Kozai mechanism where ec- 
centricity growth is possible for long-lived triple systems. If a suf- 
ficiently large value is reached, the semi-major axis and eccentric- 
ity starts decaying according to well-known relativistic expressions 
(Peters 1964). The subsequent evolution can take several forms, de- 
pending on the influence of any perturbers. 

Figure 5 (middle curve) shows the decreasing semi-major 
axis due to the full post-Newtonian treatment until coalescence at 
a coa i ~ 1.5 x 10 -11 or 6 x 10~ 4 Rq. Note the short time inter- 
val corresponding to ~ 4 x 10 3 yr required for the major part of 
this inspiral. Also plotted are the orbit-averaged solutions from Pe- 
ters (1964). Because the initial eccentricity is only available to four 

6 The requirement for parabolic escape is v csc ~ 2 km s _1 at 2r t ;<j e . 
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Figure 5. Semi-major axis as function of TV-body time. Middle curve: 
actual plot from Model BH12. Upper and lower curve: orbit-averaged so- 
lution from Peters (1964) for two initial values of eccentricity, 0.9996 and 
0.99965, where only the four first significant figures are available. The plot 
is truncated at 3 X 10 — 10 while the calculation extends to 1.5 X 10 — 11 . 
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Figure 6. Maximum eccentricity as function of time, Model BH12. Note 
the late increase of eccentricity due to Kozai cycles. 



significant figures, we have included a trial solution for a possible 
maximum value (lower curve) which brackets the actual solution. 
In view of the uncertainty, the qualitative agreement is satisfactory. 
The late stages of shrinkage are characterized by negligible dynam- 
ical perturbations. Nevertheless, in the absence of special proce- 
dures, a large number of orbits needs to be integrated. 

The corresponding plot is shown in Fig. 6 of the maximum ec- 
centricity (Heggie 1995, personal communication). This has been 
converted to mainly relativistic forrrQ and illustrates the role of the 
Kozai mechanism. In fact, a large value exceeding 0.998 was al- 
ready reached some 0.4 Myr earlier but the most significant shrink- 
age took place after the final maximum as well as the actual eccen- 
tricity exceeded 0.999. Since the eccentricity did not show any evi- 
dence of decrease during the interval, this behaviour is an example 
of so-called Kozai migration (Wu & Murray 2003). However, the 
earlier eccentricity maximum would be sufficient to ensure gravi- 
tational coalescence on a slightly longer time-scale, provided the 
favourable configuration is preserved. It can also be seen that the 



With two terms of e 2 replaced by relativistic values. 



Table 1. Summary of coalescence events. The coalescence time in Myr is 
given in Column 2, followed by the corresponding semi-major axis, eccen- 
tricity and combined mass in Mq . 
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Figure 7. Final velocity dispersion profiles as a function of central dis- 
tance, Model BH7 (upper curve) and BH8 (lower curve). The current den- 
sity centre is used as reference point and the respective times are 750 Myr 
and 1 . 1 Gyr. BHs and neutron stars have been excluded but the small mem- 
bership hardly affects the outcome. 

Einstein shift given by equation JlOt would still be small during 
this episode, consistent with the derived value TKozai — 1. 

Including all the PN terms up to 3PN demonstrates that the 
coalescence condition of equation i ll It can be reached at the maxi- 
mum eccentricity in agreement with estimates of the time-scale ©. 
This is by no means assured since it is known that the relativistic 
precession acts to de-tune eccentricity growth. In this connection, a 
discussion of the post-Newtonian modification of the Hamiltonian 
is of interest (Miller & Hamilton 2002). It turns out that the mod- 
ified expression for the maximum eccentricity still admits values 
above 0.999 for the parameters of interest here (Aarseth 2007). 

A summary of all the coalescence events is given in Table 1 . A 
total of 14 models have been investigated using the PN formulation. 
However, a few were terminated prematurely because of technical 
difficulties. Most of the other models were continued further, unless 
the total BH mass was small. Although all the cases listed were 
accepted as coalescence, some did not reach the actual end state 
of equation i ll It because one of the secondary criteria discussed 
above was satisfied. For example, in the case of Model BH1 1, the 
eccentricity was extremely large for a short time interval. Likewise, 
for Model BH14 the high eccentricity and short time-scale would 
ensure coalescence. 

There have been many observational efforts to determine the 
velocity dispersion as a function of radius, in particular related 
to the quest for discovering an intermediate mass BH in a globu- 
lar cluster. The velocity profiles of all luminous members in two 
models are shown in Fig. 7 for illustration. Here the more evolved 
Model BH8 exhibits a larger central velocity but both models are 
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Figure 8. Final density profiles, Model BH7 (upper curve) at 750 Myr and 
BH8 (lower curve) at 1 . 1 Gyr. A small number of BHs and neutron stars are 
included. 
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Figure 9. Average stellar mass in solar units as a function of central dis- 
tance. The plots for Model BH7 at 750 Myr and BH8 at 1.1 Gyr are very 
similar, with the upper stellar mass limits near 3 Mq . 



characterized by a slight velocity decrease in the innermost region. 
This deficiency may be expected to disappear as the cluster is de- 
pleted of BHs. We note that in an investigation involving a popula- 
tion of BHs, the central velocity profile (line-of-sight) showed a flat 
slope inside 1 pc (Mackey et al. 2008, Fig. 12). However, a direct 
comparison is difficult because here some 200 BHs were retained. 

The density is also of considerable interest. The final densities 
of the two models shown in Fig. 8 are quite well matched. Consis- 
tent with the small velocity decrease towards the centre, the central 
density profile has a flatter slope. Again this feature will tend to 
disappear as the post-BH core formation proceeds. 

The process of mass segregation plays an important role in the 
present work. In Fig. 9 we illustrate the average mass as a function 
of central distance in two typical models where BHs and neutron 
stars have been excluded. At the end states the maximum mass has 
been reduced to 3 Mq due to stellar evolution while the minimum 
mass used is 0.1 Mq. Although the final times differ somewhat 
(750 Myr vs 1 . 1 Gyr), the two plots are essentially indistinguishable 
and span a significant mass range. 



5 DISCUSSION 

The main purpose of this paper is to demonstrate the numerical 
capability of the code NBODY7 for studying strong interactions in 
cluster cores. It turns out that the algorithmic regularization method 
is also able to treat post-Newtonian terms with high accuracy. In 
fact, energy conservation is now substantially better than in the 
standard NBODY6 code because the dominant errors are associ- 
ated with the most strongly bound members and the ARC method 
is more accurate. This improved treatment comes at a price but al- 
ready it is feasible to investigate small globular clusters using real- 
istic parameters. 

The ARC software usually maintains high accuracy as a re- 
sult of specifying a small tolerance in the Bulirsch-Stoer integra- 
tor combined with a more efficient time transformation. A further 
advantage compared to standard NBODY6 integration with chain 
regularization is that the subsystems are treated for much longer 
times so there are less errors due to switching. In any case, a di- 
rect comparison with the chain code is only possible in the absence 
of PN terms. Moreover, the latest version of ARC admits solutions 
with only two members when PN terms are present. It is also help- 
ful that the smaller density in the inner region reduces the energy 
errors associated with the single particles. 

The onset of conditions for gravitational radiation is usu- 
ally initiated by the Kozai mechanism. Following the subsequent 
shrinkage, the inner binary orbit may become detached (or isolated) 
so that the final stages are often unperturbed dynamically. This be- 
haviour would enable a sequential analytical continuation of this 
time-consuming process, taking into account the Einstein shift of 
equation d 1 0| > or the addition of its second-order equivalent 

Stellar systems containing a population of massive objects ex- 
perience mass segregation on a short time-scale. Although the ve- 
locity kick at formation increases the velocity dispersion of neu- 
tron stars and BHs, the remaining BHs soon begin to concentrate 
towards the centre. However, the emergence of a dominant binary 
acts to deplete the central density and prevents the formation of 
other binaries, an aspect which has not been considered so far. 

The problem of black hole dynamics has been studied by many 
authors. Most studies, however, have been concerned with choices 
of the IMF which are more applicable to intermediate-mass BHs 
and hence comparisons are not appropriate. The requirement of re- 
alistic PN treatment is also rarely addressed. For these reasons, we 
only mention one recent work based on stellar mass BHs (Banerjee, 
Baumgardt & Kroupa 2010). Here an IMF truncated to low-mass 
stars (m < 1 Mq) was adopted with an addition of about 100 BHs 
of 10 Mq. Although not a fully consistent simulation, the char- 
acteristic behaviour of mass segregation was seen, together with 
ejection of massive members. Moreover, employment of the rela- 
tivistic orbit- averaged changes in semi-major axis and eccentricity 
(Peters 1964) gave rise to several cases of coalescence. 

Contrary to the evidence for BH binary formation, a full sim- 
ulation by NB0DY7 with all 140 BHs retained showed that few bi- 
naries are formed dynamically during the first Gyr when half the 
BHs have escaped. This behaviour can be understood in terms of 
an emerging binary being exposed to disruptive encounters with the 
dominant binary before becoming hard. There remains the possi- 
bility, however, that a dominant binary may be ejected temporarily 
from the core by recoil, giving the more weakly binary an opportu- 
nity of hardening. 

The question of BH coalescence by the Kozai mechanism in 

8 The analytical continuation procedure has been implemented recently. 
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long-lived triple systems has also been discussed elsewhere (Miller 
& Hamilton 2002). However, this work assumes a certain fraction 
of BH binaries in order for binary-binary collisions to take place 
and is therefore not consistent with the current findings. Given the 
relatively small number of BHs expected from a realistic IMF, the 
potential for massive binary formation is fairly limited. Moreover, 
the possible depletion due to natal velocity kicks should also be 
taken into account. 

For completeness, a runaway scenario for compact objects was 
also studied with a full post-Newtonian formulation (Kupi, Amaro- 
Seoane & Spurzem 2006). In this early PN application with the 
code NBODY6++, the treatment of all close two-body encounters 
was modified by including the effect of relativistic perturbations. 
Given the large central velocity dispersion of 4300 km s -1 , a mas- 
sive object was formed by runaway accretion after core collapse. 
We note that the Kustaanheimo-Stiefel (1965) Hermite implemen- 
tation also requires the first time derivative of the relevant PN force 
terms while in algorithmic or chain regularization only the basic 
expressions of equation I0 are needed. 



6 CONCLUSIONS 

This is a first TV-body implementation of the new algorithmic reg- 
ularization scheme (Hellstrom & Mikkola 2010). The application 
to a system containing black holes poses certain numerical prob- 
lems, mainly in connection with membership selection. Although 
the post-Newtonian terms require time-consuming calculations, it 
has proved possible to describe the full binary evolution up to co- 
alescence using realistic parameters. Further progress can be made 
taking advantage of analytical continuation of unperturbed bina- 
ries in the relativistic regime after allowing for the precession. The 
question of how to initialize the spin must also be addressed. Its 
contribution to the post-Newtonian terms enters at the order c -3 
and is therefore more important than the 2PN termjfl On the ob- 
servational side, we emphasize the presence of additional escaping 
stars of high velocity, especially in the later stages when massive 
stellar binaries are no longer present. 

The present results indicate that there is no evidence for run- 
away evolution of BH masses. Given the relatively small number of 
BHs involved, just one binary may dominate the central region. Oc- 
casionally strong interactions lead to fast ejection of other BHs and 
even the binary may escape in a slingshot event. These simulations 
show that even a small number of BHs may play an important role 
in the first Gyr of globular cluster evolution. Since the simulations 
only cover the first Gyr so far, later stages may repopulate the core 
in accordance with expectations. In the absence of a definite model 
for natal BH velocity kicks, further explorations of larger popula- 
tions are desirable. We have seen that the Kozai mechanism can be 
very effective in leading to the relativistic coalescence of BH bi- 
naries. As far as stellar mass BHs are concerned, sufficiently high 
eccentricities can be reached without the growth being de-tuned by 
precession effects. Finally, the presence of even a few stellar mass 
black holes with luminous companions may provide a new chal- 
lenge for future observations. 



7 ACKNOWLEDGMENTS 

Vital contributions to this work were made by Keigo Nitadori who 
implemented the GPU system and Seppo Mikkola who provided 
the software for algorithmic regularization. I am also indebted to 
the referee for many constructive suggestions which helped to im- 
prove the manuscript. 



REFERENCES 

Aarseth S. J„ 1971, Ap&SS, 13, 324 

Aarseth S. J., 2003a, Gravitational N-Body Simulations, Cam- 
bridge Univ. Press, Cambridge 
Aarseth S. J., 2003b, Ap&SS, 285, 367 
Aarseth S. J., 2007, MNRAS, 378, 285 
Aarseth S. J., Zare K„ 1974, Celes. Mech., 10, 185 
Banerjee S., Baumgardt H., Kroupa P., 2010, MNRAS, 402, 371 
Berentzen I., Preto M., Berczik P., Merritt D., Spurzem R., 2009, 
ApJ, 695, 455 

Blanchet L., Iyer B., 2003, Class. Quant. Grav., 20, 755 
Brockamp M., Baumgardt H., Kroupa P., 2011, MNRAS, 418, 
1308 

BulirschR., Stoer J., 1966, Num. Math., 8, 1 

Einstein A., 1915, Sitzungsber. Preuss. Acad. Wiss., 47, 2, 831 

Eldridge J., Tout C. A., 2004, MNRAS, 353, 87E 

Gieles M., Baumgardt H., Heggie D. C, Lamers J. G. L. M., 2010, 

MNRAS, 408, 16 
Harfst S, Gualandris A., Merritt D., Mikkola S., 2008, MNRAS, 

389,2 

Hellstrom C, Mikkola S., 2010, Celes. Mech. Dyn. Ast., 106, 143 
Hurley J. R., Pols O., Tout C. A., 2000, MNRAS, 315, 543 
Hurley J. R., Pols O., Tout C. A., 2002, MNRAS, 329, 897 
Iwasawa M., Funato Y., Makino J., 2006, ApJ, 651, 1059 
Kiseleva L. G., Eggleton P. P., Mikkola S., 1998, MNRAS, 300, 
292 

Kozai Y., 1962, AJ, 67, 591 
Kozai Y, 2004, Proc. Jpn. Acad., Ser. B, 80, 157 
Kroupa P., Tout C, Gilmore G, 1993, MNRAS, 262, 545 
Kupi G, Amaro-Seoane P., Spurzem R., 2006, MNRAS, 371, L45 
Kustaanheimo P., Stiefel E., 1965, J. Reine Angew. Math. 218, 
204 

Mackey A. D., Wilkinson M. L., Davies M. B., Gilmore G. F., 

2008, MNRAS, 386, 65 
Mikkola S., Aarseth S. J., 1993, Celes. Mech. Dyn. Ast., 57, 439 
Mikkola S., Aarseth S. J., 2002, Celes. Mech. Dyn. Ast., 84, 343 
Mikkola S., Merritt D., 2006, MNRAS, 372, 219 
Mikkola S., Merritt D., 2008, A. J. 135, 2398 
Mikkola S., Tanikawa K, 1999, MNRAS, 310, 745 
Miller M. C, Hamilton D. P., 2002, ApJ, 576, 894 
Mora T, Will C, 2004, Phys. Rev. D 69, 104021 
Nitadori K, 2009, Ph. D. thesis, University of Tokyo 
Nitadori K, Aarseth S. J., 2012, MNRAS, in press 
Peters P. C, 1964, Phys. Rev., 136, B1224 
SpitzerL., 1969, ApJ, 158, L139 

van Albada T. S., 1968, Bull. Astron. Inst. Neth., 19, 479 
Wu Y, Murray N., 2003, ApJ, 589, 605 
Zare, K, 1974, Celes. Mech., 10, 107 



Spin and a small merger recoil velocity have now been implemented. 



© 2011 RAS, MNRAS 000,[Tp?? 



